BEGIN {
    FS=","; 
	M= 20 # number of initial points
	N= 100 # size of one generation
	C= 0   # cluster id
	Missing="?"
	split("",W,"") # centroid weights
}
              {  gsub(/[ \t],""); gsub(/#.*/,"") }
NR==1         {  srand(Seed ? Seed : 1)
                for(I=1;I<=NF;I++) readHeader(I,$I)
			    next }
/^$/          { next }
              { for(I in Num) minMax(I,$I) }
length(W) < M { more() }
              { move(nearest()) }  
(NR % N) == 0 { less() }
END           { for(C in Last) dump(C) }

function more( i) {
	while(length(W) < M) {
		C++
		W[C]=1
    	for(i=1;i<=NR;i++) Centroid[C,i]=$i
		next
}}
function less(   sum,c,doomed,i,enough) {
	enough = rand()
	for(c in W) sum += W[c]
	for(c in W) if (W[c]/sum < enough) doomed[c]
	for(c in doomed) {
		delete W[c]
		for(i=1;i<=NR;i++) delete Centroid[c,i]
	}	
	delete Last
	for(c in W) Last[c] = W[c] = 1
}
function move(c) {
	for(i=1;i<=NF;i++) {
		if (Centroid[c,i] == Missing)
			Centroid[c,i] == $i
		else if ($i != Missing)
			Centroid[c,i] = (i in Num) ? moveNum(c,i,$i) : moveSym(c,i,$i)
	}
	W[c]++
}
function moveNum(c,col,x) {
	return (Centroid[c,col]*W[c] + x)/(W[c]+1)	
}
function moveSym(c,col,x) {
	return (rand() < 1/W[c]) ? x : Centroid[c,col]
}
function nearest(   min,c,tmp,out) {
	min = 100000
	for(c in W) 
		if((tmp = dist(c)) < min) {
			min = tmp
			out = c }
	return out
}
function dist(c,    inc,y,n, i,out) {
	for(i in Indep) {
		n++
		y = Centroid[c,i]
		inc =(i in Num) ? distNum(norm(i,$i),norm(i,y)) : distSep($i,y) 
	    out += inc^2
	}
	return sqrt(out)/sqrt(n+0.0000001)
}	
function distSym(a,b) {
	if (a == Missing || b == Missing) return 1
	return a == b ? 0 : 1
}
function distNum(a,b) {
	if (a == Missing && b == Missing) return 1
    if (a == Missing) return b < 0.5 ? 1 - b : b
	if (b == Missing) return a < 0.5 ? 1 - a : a
	return a - b
}
function minMax(col,x) {
    if(x == Missing) return 1
	if(col in Min) {if(x < Min[col]) {Min[col]=x}} else {Min[col]=x}
	if(col in Max) {if(x > Max[col]) {Max[col]=x}} else {Max[col]=x}
}
function norm(col,x) {
    if (x==Missing) return Missing
	return (x - Min[col]) / (Max[col] - Min[col] + 0.00001)
}
function readHeader(col,x,       ignore,num,dep) {
	ignore    = gsub(/\?/,"",x) 
	num       = gsub(/\$/,"",x)
	dep       = gsub(/\!/,"",x) 
	Name[col] = x
    Index[x]  = col
	if (!ignore)  {
		num ? Num[col] : Sym[col]
		dep ? Dep[col] : Indep[col]
}}
function o(a,str,     i,com) {
	com = "sort -k 2 #" rand()
	for(i in a) print str " [" i "]= "a[i] | com
	close(com)
}
function dump(c,   i) {
	print "\n Cluster #" c	
	for(i=1;i<=NF;i++)
		print "\t" Centroid[c,i]
}
